Smartphone keyboard dynamics predict affect in suicidal ideation

While digital phenotyping provides opportunities for unobtrusive, real-time mental health assessments, the integration of its modalities is not trivial due to high dimensionalities and discrepancies in sampling frequencies. We provide an integrated pipeline that solves these issues by transforming all modalities to the same time unit, applying temporal independent component analysis (ICA) to high-dimensional modalities, and fusing the modalities with linear mixed-effects models. We applied our approach to integrate high-quality, daily self-report data with BiAffect keyboard dynamics derived from a clinical suicidality sample of mental health outpatients. Applying the ICA to the self-report data (104 participants, 5712 days of data) revealed components related to well-being, anhedonia, and irritability and social dysfunction. Mixed-effects models (55 participants, 1794 days) showed that less phone movement while typing was associated with more anhedonia (β = −0.12, p = 0.00030). We consider this method to be widely applicable to dense, longitudinal digital phenotyping data.


Model order analysis
In the main text, we discussed the most parsimonious 5-component solution.Different ICA model orders might be just as valid, however.To examine the behaviour of the independent components across multiple model orders, we reran our analyses with ten and twenty components and correlated the independent components of those solutions with the independent components of the 5-component solution.

Ten components:
The mixing matrix of a 10-component ICA solution is given in Supplementary Figure 1 and the cross-correlations between the independent components of this solution and the original 5-component solution are given in Supplementary Figure 2. To facilitate the comparison between the ICs of these solutions, here we will refer to the components by the prefix ICX, where the subscript X indicates the model order (e.g., IC5 for an IC of the 5-component solution).IC5 1 correlates most strongly with IC10 4, which is expected given the large loadings on the well-being variables in their respective mixing matrices.(Note that the sign of the correlations does not matter, as the independent components themselves are only defined up to a multiplicative sign. 1 ) IC5 2, on the other hand, shows strongest correlations with both IC10 2 and IC10 5.The correlation with IC10 2 can be explained by its large loadings on the anhedonia items, and the correlation with IC10 5 can be explained by a comparable loading polarity pattern (i.e., negative loadings on positive variables, positive loadings on negative variables, and relatively large loadings on the DRSP variables).IC5 3 shows the strongest correlation with IC10 4, but also a moderate correlation with IC10 8. We see the same pattern as with IC5 2: IC10 8 has loadings for which the size  1.
corresponds to the loadings of IC5 3, while IC10 4 displays a loading polarity pattern similar to that of IC5 3. It is not unlikely that IC5 3 has split up into these two IC10 components.IC5 4 shows a strong correlation with IC10 6, which is not surprising given their similar loading values.IC5 5, finally, shows its strongest correlation with IC10 1. Their loadings are somewhat comparable, with positive loadings on the agitation items (BAM), negative loadings on the irritability items (BITe), and mixed loadings on the DRSP items.
As for the mixed-effects models of the 10-component solution, none of the effects survived Bonferroni correction (see Supplementary Table 2).We point out, however, that the nominally significant, negative effect of phone movement on IC10 5 is consistent with what we found for the anhedonia component of the 5-component solution.
All models except the one for IC10 8 showed no signs of heteroskedasticity.Most models, except those for IC10 5 and 10, showed slight departures from normal residuals.Participant-level random effects of the models for IC10 1 and 6-10 displayed small to medium deviations from normality.Week-within-participant-level random effects showed small deviations from normality for all models except for the one for IC10 10.

Twenty components:
The mixing matrix of the 20-component ICA solution is given in Supplementary Figure 3 and the cross-correlations between its components and the components of the 5-component solution are given in Supplementary Figure 4. IC5 1 and IC5 2 correlate most strongly with IC20 7. The 20-component mixing matrix shows that IC20 7 has relatively strong loadings on both the well-being and anhedonia items, suggesting that it is a recombination of IC5 1 and IC5 2. Our phone movement sensitivity analysis (described below) shows results consistent with this suggestion.IC5 3 shows its strongest correlations with IC20 10 (similar loading pattern) and IC20 13 (similar loading magnitude), which is the same behaviour we found for the 10-component solution.IC5 4 has a strong correlation with IC20 2, which is (again) unsurprising because they display very similar loadings.IC5 5 displays its strongest correlation for IC20 15: Both components again show positive agitation loadings, (mostly) negative irritability loadings, and mixed DRSP loadings.
Consistent with what we found in our sensitivity analysis, the mixed-effects model of the general affect component of the 20-component solution (IC20 7) shows a significant association with phone movement after Bonferroni correction (β = -0.14, p = 0.00011).The sign of this effect is also consistent with the effect of phone movement on the anhedonia component in the 5-component solution (IC5 2): More movement predicts lower ratings on LackingInterest, Anhedonia, and Unmotivated.1.
Only the model of IC20 13 showed signs of heteroskedasticity, but several other models (IC20 1, 6, 18, and 20) displayed structure in their residuals versus fitted-values plots.All models except the ones for IC20 3, 7, and 18 showed small to medium deviations from residual normality.Participant-level random effects deviated moderately from normality for all models except the ones for IC20 9, 10, 11, and 18. Week-within-participant-level random effects deviated slightly from normality for all models except the ones for IC20 7, 18, and 20.

Phone movement sensitivity analysis
Because the FastICA algorithm starts with a random estimate, its final solution differs from run to run. 1 We conducted a sensitivity analysis to estimate the impact of these differences on the effect of phone movement.More specifically, we ran the fastICA function 100 times on our data, manually examined all ICA solutions to select the component that best matched the anhedonia component and constructed a mixed-effects model for every selected component.This procedure resulted in 100 effect estimates of phone movement.
The ICA solutions displayed a dichotomy; an example from each class of solutions is shown in Supplementary Figure 5. 69% of the solutions featured a general affect component, with large loadings on the positive affect variables and large, opposite loadings on the negative affect variables, especially for the DRSP variables.The remaining 31% showed the well-being and anhedonia components previously encountered in the main text, which seem a split of the general affect component.
In the ICA solutions where the anhedonia component did not appear, we created a model of the general component instead.All 100 models reported significant associations with phone movement after Bonferroni correction within a single sensitivity iteration (-0.14 £ β £ 0.14, all p £ 0.00031).The sign of the β values were flipped depending on whether the sign of their associated independent component was also flipped.The size of the β and p values depended on the class of the ICA solution: If the ICA solution contained the general component, |β| = 0.14 and p < 0.0001.Otherwise, |β| = 0.12 and p £ 0.00031.These results can be interpreted as a dilution of the effect when the general component is split into a well-being and anhedonia one.
Even though the general affect solution was more prevalent in our sensitivity analysis, we still decided to discuss the solution with separate well-being and anhedonia components in our main text because 1) the well-being plus anhedonia solution appeared first in our analyses and 2) the splitting of general affect into well-being and anhedonia provides a more fine-grained view of what drives the association with phone movement.).IC = independent component.For the questionnaire abbreviations, please refer to Supplementary Table 1.

Contiguous data analysis
In the main text, we ran the ICA on all available self-report data.This meant that we also included data from periods with high proportions of missing data, due to which the data stream is not contiguous but fragmented.
Temporal ICA is well suited to handle fragmented data since optimising for statistical independence requires assessing the probability density of a source process.ICA does assume, however, that the data were generated from a stationary process.This may not be the case when the distribution underlying data from the periods of fragmentation are very different from the distribution of the contiguous data.Participants might, for example, be less inclined to fill out the self-report items when they are having a bad day and bias the self-report surveys to only measure good days, or vice versa.This non-stationarity could, in turn, affect our ICA results.To test this hypothesis, we reran our sensitivity analysis with the subset of the data that conformed to a contiguity constraint: Self-report data had to be present for at least seven contiguous days.Blocks smaller than seven days would be discarded.
Supplementary Figure 6 shows what the data inclusion patterns look like when the contiguity constraint is applied.
Because some of the self-report data is excluded, some of the BiAffect data also had to be excluded due to the complete-case requirements of the linear regression models.In total, 4215 days' worth of self-report data (98 In the left panel, all BiAffect data that is missing or falls outside the included range of self-report data is marked as excluded. participants) were fed into the ICA, and 1454 days' worth of BiAffect and ICA component data (47 participants) were fed into the mixed-effects models.
To determine the behaviour of the ICA with contiguous data, we ran the same sensitivity analysis as for the original data.The contiguous ICA solutions showed the same dichotomy as the fragmented ones, but in different proportions: 87% (previously 69%) of the ICA solutions showed the general affect component, while the remaining 13% (previously 31%) showed the split into the well-being and anhedonia components.The solutions with the general component showed a significant effect of phone movement after Bonferroni correction (after rounding, all |β| = 0.12, all p £ 0.00020); the solutions with the anhedonia component did not (|β| £ 0.070, 0.41 £ p £ 0.53).The fact that the proportion of solutions with an anhedonia component has decreased could be interpreted as the consequence of some non-stationarity caused by the fragmented data.Apparently, the fragmented data provide more evidence for separate well-being and anhedonia components, perhaps because the anhedonia items featured more prominently in the fragmented periods.It is, however, also possible that we found the two components more often in the fragmented case simply because that case has more data available to discriminate the two components.The non-significant effects of phone movement on the anhedonia component in the contiguous case contrast with the significant effects in the fragmented case.This contrast can partially be explained by a reduction in statistical power due to the lowered number of data points.
As we pointed out above, if the distributions underlying our data were indeed non-stationary and the periods of missingness in our self-report data were non-random, our entire analysis could be biased towards the days in which participants felt relatively good or bad.Unfortunately, solving the problem of data that are not missing at random (MNAR) is only possible when collecting additional data or making assumptions about the missing data mechanisms. 2 We currently have too little knowledge about what happens to our participants during the periods of missingness and what effect that would have on both the self-report and keyboard dynamics data to properly impute the missing data.For that reason, we decided to do a complete case analysis instead, accepting the potential bias that would introduce.The model order analysis was also repeated for the contiguous case.For brevity, we do not present the results here, but the patterns in the higher-order ICA solutions are similar to those found in the fragmented case.

Within-participants mean-centring
In the main text, we found a component with moderate to low negative loadings on all self-report items (IC 3).To examine whether this component might represent a mean offset of self-report ratings for some of our participants, we averaged the self-report ratings per participant and compared them to the component values.More specifically, we first averaged self-report item ratings within participant and within item, and then averaged again within participants but across items.IC 3 values were also averaged within participants.The averaged self-report and

Forwards-fitting of random effects
We tested the benefit of random slopes by running a forwards-fitting procedure on the 5-component mixed-effects models as specified in the main text.This procedure would add one random slope variance parameter at a time, either on level 1 (participant) or level 2 (week within participant), constraining the random effects covariance matrix to be diagonal.Likelihood ratio tests were conducted for all pairs of nested models to decide whether the additional slope parameter improved the model significantly.
The resulting models are shown in Supplementary Table 5.Notice that for the models for IC 2 and 3, the procedure did not converge to a single largest model.For IC 5, there were no significant improvements over the base model.None of the models had a large impact on the fixed-effects estimates, which we do not show for brevity.In other words, the conclusions we drew from the models in the main text remained unchanged.

2 :
Cross-correlations between the independent component time series of the 5-and 10component solutions.IC = independent component.
Cross-correlations between the independent component time series of the 5-and 20component solution.IC = independent component.

Supplementary Figure 5 :
Example mixing matrices of the two classes of ICA solutions.The left panel shows a solution which features the general affect component (IC 1), while the right panel shows the anhedonia component (IC 5

Supplementary Figure 6 :
Missing data patterns with the contiguity constraint.In the right panel, all self-report data that are missing or do not comply with the contiguity constraint are marked as excluded (red strikethrough).
Supplementary Figure 7: A participant's average daily rating and IC 3 value over time.

Supplementary Figure 3: Mixing matrix of the 20-component solution
IC 11 IC 12 IC 13 IC 14 IC 15 IC 16 IC 17 IC 18 IC 19 IC 20 .IC = independent component.For the questionnaire abbreviations, please refer to Supplementary Table

Supplementary Table 2: Model estimates for the 10-component ICA solution.
Every IC corresponds to a separate model.p' indicates uncorrected p values.IC = independent component; IKD = inter-key delay; MAD = mean absolute deviation.(Continues on next page.)

Table 3 : Model estimates for the 20-component ICA solution.
Every IC corresponds to a separate model.p' indicates uncorrected p values.IC = independent component; IKD = inter-key delay; MAD = mean absolute deviation.(Continues on next page.)

Table 5 Random slopes that were significant additions to their nested models
. IC = independent component; IKD = inter-key delay.